Surowy zbiór danych był na tyle duży, że problematyczne było załadowanie go w całości do pamięci RAM. Dopiero usunięcie pewnej części zbędnych danych jeszcze przed wczytaniem pozwoliło na wczytanie go do data frame’a. Mimo usunięcia sporej liczby krotek i atrybutów podczas wstępnego przetwarzania, nadal w fazie uczenia maszynowego mieliśmy doczynienia z ponad 350000 obserwacjami, po 383 atrybuty każda. Było to zbyt dużo dla większości algorytmów, aby być w stanie zbudować model korzystając z pozostałej wolnej pamięci RAM. W tym wypadku pierwotnie zdecydowano się na podejście przyrostowe i budowanie modelu stopniowo z na tyle małych wycinków zbioru danych, aby było możliwe przeprowadzenie uczenia. Później jednak wybrano selekcję atrybutów jak sposób na rozwiązanie problemów z pamięcią. Przy okazji przyspieszyło to znacznie same obliczenia.
library(dplyr)
library(ggplot2)
library(caret)
library(DT)
library(reshape2)
library(tidyr)
library(plotly)
library(biglm)
library(xgboost)set.seed(95)Usunięcie atrybutów redundantnych, pustych, o zbyt złożonej strukturze do analizy. Decyzję podjęto na podstawie wstępnego przejrzenia pliku all_summary.csv zawierającego cały zbiór danych. Pozwala to zmniejszyć rozmiar wczytywanego pliku, a w konsekwencji rozmiar wynikowego data frame’a. Dzięki temu zabiegowi jest on w stanie zmieścić się w pamięci RAM.
if (!("data.csv" %in% dir()))
system("cut -d';' -f1-3,27,399 --complement all_summary/all_summary.csv > data.csv")Wczytanie danych i usunięcie kolumn które nie wezmą udziału w dalszym przetwarzaniu.
df <- read.table("data.csv", nrows = 75000, header = T, sep=";", comment.char = "")
classes <- sapply(df, class)
df <- read.table("data.csv", nrows = 591043, header = T, sep=";", comment.char = "", colClasses = classes)
colToRemove <- c("local_BAa", "local_NPa", "local_Ra", "local_RGa", "local_SRGa", "local_CCSa", "local_CCPa", "local_ZOa", "local_ZDa", "local_ZD_minus_a", "local_ZD_plus_a", "local_res_atom_count", "local_res_atom_non_h_occupancy_sum", "local_res_atom_non_h_electron_occupancy_sum", "local_res_atom_C_count", "local_res_atom_N_count", "local_res_atom_O_count", "local_res_atom_S_count", "dict_atom_C_count", "dict_atom_N_count", "dict_atom_O_count", "dict_atom_S_count", "fo_col", "fc_col", "weight_col", "grid_space", "solvent_radius", "solvent_opening_radius", "part_step_FoFc_std_min", "part_step_FoFc_std_max", "part_step_FoFc_std_step")
df <- df[, !(names(df) %in% colToRemove)]Usunięcie ze zbioru danych ktrotek, które posiadają wskazane wartości atrybutu res_name
df <- filter(df, !res_name %in% c("UNK", "UNX", "UNL", "DUM", "N", "BLOB", "ALA", "ARG", "ASN", "ASP", "CYS", "GLN", "GLU", "GLY", "HIS", "ILE", "LEU", "LYS", "MET", "MSE", "PHE", "PRO", "SEC", "SER", "THR", "TRP", "TYR", "VAL", "DA", "DG", "DT", "DC", "DU", "A", "G", "T", "C", "U", "HOH", "H20", "WAT"))Na tym etapie przetwarzania w zbiorze danych występuje 40036 krotrek o conajmniej jednej wartości pustej. Przekłada się to na nieco mniej niż 7% całego zbioru. Ponieważ dalej chcemy rozważać tylko 50 najczęstszych wartości atrybutu res_name sprawdzimy, czy usunięcie wszytskich krotek z wartościami pustymi wpłynie na ten wybór. Okazuje się, że po zostawieniu jedynie kompletnych krotek nadal wybierzemy te same 50 wartości res_name. W związku z tym decydujemy się na tę obsługę wartości pustych.
sumna <- sapply(df, function(x) sum(is.na(x)))
sumna <- sumna[sumna > 0]
print(max(sumna))## [1] 40036print(max(sumna)/length(df[[1]]))## [1] 0.06839797dfComplete <- df[complete.cases(df),]
countn <- (df %>% group_by(res_name) %>% summarise(n()) %>% rename(sum = "n()") %>% filter(!is.na(res_name)) %>% arrange(desc(sum)))[1:50,1]
countn2 <- (dfComplete %>% group_by(res_name) %>% summarise(n()) %>% rename(sum = "n()") %>% arrange(desc(sum)))[1:50,1]
length(intersect(countn$res_name, countn2$res_name))## [1] 50df <- dfComplete
rm(dfComplete)Po wstępnym przetworzeniu zbiór danych składa się z 525666 krotek, każda zawiera 388 atrybutów. 4 z nich są typu factor a pozostałe 384 to liczby, gdzie 51 to liczby całkowite a 333 zmiennoprzecinkowe.
print(length(df[[1]]))## [1] 525666print(length(df))## [1] 388knitr::kable(table(sapply(df, class)))| Var1 | Freq |
|---|---|
| factor | 4 |
| integer | 51 |
| numeric | 333 |
Sprawdzenie liczności każdej z 50 klas.
df <- filter(df, res_name %in% countn2$res_name)
countn <- df %>% group_by(res_name) %>% summarise(n()) %>% rename(sum = "n()") %>% arrange(desc(sum))
df$res_name <- factor(df$res_name)
prettyTable(countn)Po wyliczeniu korelacji między atrybutami liczbowymi okazuje się, że wiele z nich jest bardzo silnie skorelowanych.
correl <- cor(df[, sapply(df, is.numeric)])
correl[lower.tri(correl)] <- NA
correl <- melt(correl)
correl <- correl[complete.cases(correl), ] %>% filter(Var1 != Var2) %>% mutate(absValue = abs(value)) %>% arrange(desc(absValue))
prettyTable(correl)cantUse <- c("dict_atom_non_h_count", "dict_atom_non_h_electron_sum", "local_res_atom_non_h_electron_sum", "local_res_atom_non_h_count")
colForElectron <- (correl %>% filter(Var1 == "local_res_atom_non_h_count", abs(value) > 0.6))$Var2
colForElectron <- as.vector(colForElectron[!colForElectron %in% cantUse])
colForAtom <- (correl %>% filter(Var1 == "local_res_atom_non_h_electron_sum", abs(value) > 0.6))$Var2
colForAtom <- as.vector(colForAtom[!colForAtom %in% cantUse])
rm(correl)Zarówno rozkład liczby atomów jak i elektronów skupia się głównie na mniejszych wartościach, stosunkowo niewielkich w porównaniu do zakresu wartości na tych atrybutach.
ggplotly(ggplot(df, aes(local_res_atom_non_h_count)) + geom_density())ggplotly(ggplot(df, aes(local_res_atom_non_h_electron_sum)) + geom_density())Ze względu na stosunkowo duży zakres wartości obu atrybutów, zdecydowano się bazować na błędzie względnym zamiast bezwzględnego. Podejście te ma tą własność, że pomyłka o 2 elektrony będzie dużym błędem gdy rzeczywista wartość to 1, natomiast znacznie mniejszym przy rzeczywistej wartości 400.
knitr::kable((df %>% group_by(res_name) %>% summarise(mean(abs(1-local_res_atom_non_h_count/dict_atom_non_h_count))) %>% rename(sum = 2) %>% arrange(desc(sum)))[1:10,] %>% rename("średni błąd względny liczby atomów" = 2))| res_name | średni błąd względny liczby atomów |
|---|---|
| 1PE | 0.1671000 |
| MLY | 0.1093553 |
| SEP | 0.0909731 |
| PG4 | 0.0799353 |
| MAN | 0.0723490 |
| NAG | 0.0653264 |
| CLA | 0.0550504 |
| PLP | 0.0539990 |
| PGE | 0.0402994 |
| COA | 0.0392920 |
knitr::kable((df %>% group_by(res_name) %>% summarise(mean(abs(1-local_res_atom_non_h_electron_sum/dict_atom_non_h_electron_sum))) %>% rename(sum = 2) %>% arrange(desc(sum)))[1:10,] %>% rename("średni błąd względny liczby elektronów" = 2))| res_name | średni błąd względny liczby elektronów |
|---|---|
| 1PE | 0.1688296 |
| MLY | 0.1267699 |
| SEP | 0.0914773 |
| MAN | 0.0826846 |
| PG4 | 0.0798827 |
| NAG | 0.0760858 |
| PLP | 0.0583479 |
| CLA | 0.0527664 |
| PGE | 0.0396442 |
| COA | 0.0389528 |
W tej sekcji nastąpi próba obliczania liczby atomów i elektronów, oraz przewidzenia wartości atrybutu res_id. ##Regresja liczby atomów
moreColToRemove <- c("dict_atom_non_h_count", "dict_atom_non_h_electron_sum", "pdb_code", "res_id", "chain_id")
df <- df[, !(names(df) %in% moreColToRemove)]
forTrain <- createDataPartition(df$local_res_atom_non_h_count, p= .95, list=F)
formula <- as.formula(paste("local_res_atom_non_h_count ~ ", paste(colForAtom, collapse = " + ")))
split = 20
colForAtom <- append(colForAtom, "local_res_atom_non_h_count")
df2 <- df[, names(df) %in% colForAtom]
model <- lm(formula, data=df2[forTrain, ])
y <- predict(model, df2[-forTrain, ])## Warning in predict.lm(model, df2[-forTrain, ]): prediction from a rank-
## deficient fit may be misleadingdiff <- y-df[-forTrain, "local_res_atom_non_h_count"]
diff <- diff*diff
rmse <- sqrt(mean(diff))
r2 <- cor(y, df[-forTrain, "local_res_atom_non_h_count"])^2
print(rmse)## [1] 9.382383print(r2)## [1] 0.4902731rm(model, y, diff)forTrain <- createDataPartition(df$local_res_atom_non_h_electron_sum, p= .95, list=F)
formula <- as.formula(paste("local_res_atom_non_h_electron_sum ~ ", paste(colForElectron, collapse = " + ")))
split = 20
colForElectron <- append(colForElectron, "local_res_atom_non_h_electron_sum")
df2 <- df[, names(df) %in% colForElectron]
model <- biglm(formula, df2[forTrain, ])
y <- predict(model, df2[-forTrain, ])
diff <- y-df[-forTrain, "local_res_atom_non_h_electron_sum"]
diff <- diff*diff
rmse <- sqrt(mean(diff))
r2 <- cor(y, df[-forTrain, "local_res_atom_non_h_electron_sum"])^2
print(rmse)## [1] 65.45624print(r2)## [,1]
## [1,] 0.4620111rm(model, y, diff)Regresja nie okazała się być bardzo dokładną, jednak może to być spowodowane prostotą modelu linowego. ##Klasyfikacja W celu poprawienia jakości klasyfikacji, a przede wszystkim czasu uczenia i zajętości pamięciowej dokonana została selekcja atrybutów. Najpierw usunięto atrybuty o zerowej wariancji a następnie te, które były silnie skorelowane z innymi. Jako klasyfikatora użyto xgboost.
colToRemove <- c("local_res_atom_non_h_count", "local_res_atom_non_h_electron_sum")
df <- df[, !(names(df) %in% colToRemove)]
df <- df[,sapply(df, function(v) !is.numeric(v) || (var(v)!=0))]
correl <- cor(df[, sapply(df, is.numeric)])
correl[lower.tri(correl, diag = T)] <- 0
correl <- apply(correl, 2, max)
correl <- correl[correl > 0.8]
df <- df[, !(names(df) %in% names(correl))]
df$res_name <- factor(df$res_name)
forTrain <- createDataPartition(df$res_name, p=.7, list = F)
dtest <- xgb.DMatrix(data = as.matrix(df[-forTrain, names(df) != "res_name"]), label = df[-forTrain, "res_name"])
dtrain <- xgb.DMatrix(data = as.matrix(df[forTrain, names(df) != "res_name"]), label = df[forTrain, "res_name"])
watchlist <- list(train=dtrain, test=dtest)
bstDMatrix <- xgb.train(data = dtrain, max.depth = 10, eta = 1, nthread = 4, nrounds = 1, objective = "multi:softmax", tree_method="approx", num_class=51, watchlist=watchlist)## [23:50:39] Tree method is selected to be 'approx'
## [1] train-merror:0.359822 test-merror:0.400867